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t~ 5 Abstract. Recent works on cost based relaxations have improved Con- 

straint Programming (CP) models for the Traveling Salesman Problem 

t-H (TSP). We provide a short survey over solving asymmetric TSP with 

CP. Then, we suggest new implied propagators based on general graph 
properties. We experimentally show that such implied propagators bring 
robustness to pathological instances and highlight the fact that graph 
structure can significantly improve search heuristics behavior. Finally, 
. we show that our approach outperforms current state of the art results. 

o 

1 Introduction 

> 

Given a n node, m arc complete directed weighted graph G = (V, A, f : A — »• W), 

the Asymmetric Traveling Salesman Problem [T] (ATSP) consists in finding a 
partial subgraph G = (V, A' , /) of G which forms a Hamiltonian circuit of min- 
imum cost. This NP-hard problem is one of the most studied by the Operation 

\£) Research community. It has various practical applications such as vehicle routing 

problems of logistics, microchips production optimization or even scheduling. 
The symmetric TSP is well handled by linear programming techniques pQ. 
. . However, such methods suffer from the addition of side constraints and asymmet- 

ric cost matrix, whereas constraint programming models do not. Since the real 

^> world is not symmetric and industrial application often involve constraints such 

as time windows, precedences, loading capacities and several other constraints, 
improving the CP general model for solving the ATSP leads to make CP more 
competitive on real world routing problems. Recent improvements on cost based 
relaxations [J] had a very strong impact on the ability of CP technologies to solve 
the TSP. In this paper, we investigate how the graph structure can contribute 
to the resolution process, in order to tackle larger instances. For this purpose, 
we developed usual and original filtering algorithms using classical graph struc- 
tures, such as strongly connected components or dominators. We analyzed their 
behavior both from a quantitative (time complexity) and a qualitative (con- 
sistency level) point of view. Also, we experimentally show that such implied 
propagators bring robustness to hard instances, and highlight the fact that the 



graph structure can significantly improve the behavior of search heuristics. Our 
main contribution is both a theoretical and an experimental study which lead 
to a robust model that outperforms state of the art results in CP. 

This paper is divided into six main parts. Section[2]provides some vocabulary 
and notations. Section [3] discusses the state of the art implied constraints. Next, 
we show in Section [4] how the reduced graph can provide useful information for 
pruning and improving existing models. In Section [5] we provide some improve- 
ments about the implementation of the Held and Karp method within a directed 
context. Section [6] shows an experimental study on the ATSP and some openings 
about its symmetric version (TSP). Section [7] concludes the paper with several 
perspectives. 

2 Background 

Let us consider a directed graph G = (V,A). A Strongly Connected Component 
(SCC) is a maximal subgraph of G such that for each pair of nodes {a, b} G V 2 , 
a path exists from a to b and from b to a. A reduced graph Gr = (Vr, Ar) of a 
directed graph G represents the SCC of G. This graph is obtained by merging 
the nodes of G which are in the same SCC and removing any loop. Such a 
graph is unique and contains no circuit. We link G and Gr with two functions: 
sccOf : V -> Vr and nodesOf : Vr -> V v . The method sccOf can be represented 
by one n-size integer array. Also, since each node of V belongs to exactly one 
SCC of Vr, the method nodesOf can be represented by two integer arrays: the 
first one represents the canonical element of each SCC while the second one 
links nodes of the same SCC, behaving like a linked list. Those two arrays have 
respectively size of ur and n, where ur = | Vr]. The transitive closure of G 
is a graph Gtc = (V,Atc) representing node reachability in G, i.e. such that 
(i, j) £ Atc if and only if a path from i to j exists in G. 

In a CP context a Graph Variable can be used to model a graph. Such 
a concept has been introduced by Le Pape et al. [21] and detailed by Regin 
[28] and Dooms et al. [H]. We define a graph variable GV by two graphs: the 
graph of potential elements, Gp = (Vp,Ap), contains all the nodes and arcs 
that potentially occur in at least one solution whereas the graph of mandatory 
elements, Gm = (Vm, Am), contains all the nodes and arcs that occur in every 
solution. It has to be noticed that Gm C Gp C G. During resolution, decisions 
and filtering rules will remove nodes/arcs from Gp and add nodes/arcs to Gm 
until the Graph Variable is instantiated, i.e. when Gp — Gm- It should be 
noticed that, regarding the TSP, Vp = Vm = V, so resolution will focus on Am 
and Ap: branching strategies and propagators will remove infeasible arcs from 
Ap and add mandatory arcs of Ap into Am ■ 



3 Related Work 



This section describes the state of the art of existing approaches for solving 
ATSP with CP. We distinguish the structural filtering, which ensures that a 
solution is a Hamiltonian path, from cost based pruning, which mainly focus on 
the solution cost. Then, we study a few representative branching heuristics. 

Given, a directed weighted graph G — (V,A,f), and a function / : A — > R, 
the ATSP consists in finding a partial subgraph G' = (V, A' , /) of G which forms 
a Hamiltonian circuit of minimum cost. A simple ATSP model in CP, involving 
a graph variable GV, can basically be stated as minimizing the sum of costs 
of arcs in the domain of GV and maintaining GV to be a Hamiltonian circuit 
with a connectivity constraint and a degree constraint (one predecessor and one 
successor for each node). However, it is often more interesting to convert such 
a model in order to find a path instead of a circuit [15125) . Our motivation for 
this transformation is that it brings graph structure that is more likely to be 
exploited. 

In this paper, we consider the ATSP as the problem of finding a minimum 
cost Hamiltonian path with fixed start and end nodes in a directed weighted 
graph. In the following, s, e £ V respectively denote the start and the end of 
the expected path, s and e are supposed to be known. They can be obtained by 
duplicating any arbitrary node, but it makes more sense to duplicate the node 
representing the salesman's home. 

3.1 Structural filtering algorithms 

Our formulation of the ATSP involves the search of a path instead of a circuit, the 
degree constraints has thus to be stated as follows: For all v £ V\{e}, Sq,(v) = 1 
and for any v £ V\{s},5q,(v) = 1, where Sq,(v) (respectively Sq,(v)) denotes 
the number of successors (respectively predecessors) of v. Extremal conditions, 
being Sq, (e) = 8q, (s) = 0, are ensured by the initial domain of the graph vari- 
able. An efficient filtering can be obtained with two special purpose incremental 
propagators. One reacts on mandatory arc detections: whenever arc (u, v) is en- 
forced, other outgoing arc of u and ingoing arcs of v can be pruned. The other 
reacts on arc removals: whenever a node has only one outgoing (or ingoing) 
potential arc left, this arc is mandatory and can be enforced. A higher level of 
consistency can be achieved by using a graph-based AllDifferent constraint 
maintaining a node-successor perfect matching |27j . Deleting circuits is the sec- 
ond important structural aspect of the TSP. Caseau and Laburthe [7] suggested 
the simple and efficient NoCycle constraint to remove circuits of the graph. Their 
fast incremental algorithm is based on the subpaths fusion process. It runs in 
constant time per arc enforcing event. The conjunction of this circuit elimination 
constraint and the above degree constraints is sufficient to guarantee that the 
solution is a Hamiltonian path from s to e. 

However, other implied constraints provide additional filtering that may help 
the resolution process. For instance, Quesada [26 suggested the general propaga- 
tor DomReachability which maintains the transitive closure and the dominance 



tree of the graph variable. However, its running time, 0(nm) in the worst case, 
makes it unlikely to be profitable in practice. A faster constraint, also based on 
the concept of dominance, is the Arborescence constraint. It is nothing else but 
a simplification of the Tree constraint [5] recently improved to a 0(n + m) worst 
case time complexity |12j . Given a graph variable GV and a node s, such a con- 
straint ensures that GV is an arborescence rooted in node s. More precisely, it 
enforces GAC over the conjunction of the following properties: Gp has no circuit, 
each node is reachable from s and, each node but s has exactly one predecessor. 
Such a filtering can also be used to define the AntiArborescence by switching 
s with e and predecessors with successors. 

A dual approach consists in assigning to each node its position in the path. 
In such a case, the position of a node is represented by an integer variable 
with initial domain [0,n — 1]. Positions are different from a node to another, 
which can be ensured by an AllDifferent constraint. Since the number of 
nodes is equal to the number of positions, the bound consistency algorithm of 
AllDifferent constraint only requires 0(n) time. Plus, a channeling has to 
be done between the graph variable and position variables. Such a channeling 
requires 0(n + m) worst case time. In particular, lower bounds of positions are 
adjusted according to a single Breadth First Search (BFS) of Gp(s) while upper 
bounds of positions are shortened by a BFS of Gp 1 (e). It has to be noticed 
that this approach is related to disjunctive scheduling [33] : nodes are tasks of 
duration 1 which are executed on the same machine. The structure of the input 
graph creates implicit precedence constraints. 

Finally, some greedy procedures based on the finding of cuts have been sug- 
gested in the literature: Bcnchimol et al. enforce some cut-sets of size two [I] 
while Kaya and Hooker use graph separators for pruning JT] . The drawback of 
such methods is that they provide no level of consistency. 

3.2 Cost-based filtering algorithms 

CP models often embed relaxation based constraints, to provide inference from 
costs. Fischetti and Toth [T5] suggested a general bounding procedure for com- 
bining different relaxations of the same problem. 

The most natural relaxation is obtained by considering the cheapest outgoing 
arc of each node: LB tr iviai = J2uev\{e} mm {/( u i v )\( u j v ) € Ap}. Such a lower 
bound can be computed efficiently but it is in general relatively far from the 
optimal value. 

A stronger relaxation is the weighted version of the AllDifferent con- 
straint, corresponding to the Minimum Assignment Problem (MAP). It requires 
0(n(m + nlogn)) time [23] to compute a first minimum cost assignment but 
then 0(n 2 ) time [6J to check consistency and filter incrementally Some inter- 
esting evaluations are provided by [IB], but are mainly related to the TSP with 
time windows constraints. 

A widely exploited subproblem of the ATSP is the Minimum Spanning Tree 
(MST) problem where the degree constraint and arc direction are relaxed. We 
remark that a hamiltonian path is a spanning tree and that it is possible to 



compute a MST with a degree restriction at one node [17] , A MST can be 
computed in two ways. The first one is Kruskal's algorithm, which runs in 
0(am) worst case time, where a is the inverse Ackermann function, but re- 
quires edges to be sorted according to their weights. Sorting edges can be done 
once and for all in 0(m log m) time. The second option is to use Prim's algo- 
rithm which requires 0(m+n\ogn) time with Fibonacci heaps [17] or 0{m log n) 
time if binomial heaps are used instead. Based on Kruskal's algorithm, Regin 
et al. |29l30j made the Weighted Spanning Tree constraint which ensures con- 
sistency, provides a complete pruning and detects mandatory arcs incremen- 
tally, within 0(am) time. Dooms and Katriel [10111) presented a more complex 
Minimum Spanning Tree constraint which maintains a graph and its spanning 
tree, pruning according to King's algorithm |22j . 

An improvement of the MST relaxation is the approach of Held and Karp [TH] , 
adapted for CP by Benchimol et al. [3] . It is the Lagrangian MST relaxation with 
a policy for updating Langrangian multipliers that provides a fast convergence. 
The idea of this method is to iteratively compute MST that converge towards a 
path by adding penalties on arc costs according to degree constraints violations. 
It must be noticed that since arc costs change from one iteration to another, 
Prim's algorithm is better than Kruskal's which requires to sort edges. Moreover, 
to our knowledge neither algorithm can be applied incrementally. 

A more accurate relaxation is the Minimum Spanning Arborescence (MSA) 
relaxation, since it does not relax the orientation of the graph. This relaxation 
has been studied by |14I15| who provide a 0(n 2 ) time filtering algorithm based 
on primal/dual linear programs. The best algorithm for computing a MSA has 
been provided by Gabow et al. [17] . Their algorithm runs in 0(m + nlogn) worst 
case time, but it does not provide reduced costs that are used for pruning. Thus, 
it could be used to create a Minimum Spanning Arborescence constraint with a 
0(m + nlogn) time consistency checking but the complete filtering algorithm 
remains in 0(n 2 ) time. The Lagrangian MSA relaxation, with a MSA computa- 
tion based on Edmonds' algorithm, has been suggested in |7|. This method was 
very accurate but unfortunately unstable. Also, Benchimol et al. [3] report that 
the MSA based Held and Karp scheme lead to disappointing results. 

3.3 Branching heuristics 

Branching strategies forms a fundamental aspect of CP which can drastically re- 
duce the search space. We study here dedicated heuristics, because the literature 
is not clear about which branching heuristic should be used. 

Pesant et al. have introduced Sparse heuristic [25 which has the singularity of 
considering occurrences of successors and ignoring costs. In this way, this heuris- 
tic is based on the graph structure. It behaves as following: First, it selects the set 
of nodes X with no successor in Gm and the smallest set of successors in Gp. Sec- 
ond, it finds the node x G X which maximize J^V y - )l£Ap \{(z,y) G Ap\z e X}\. 
The heuristic then iterates on x's successors. This process is optimized by per- 
forming a dichotomic exploration of large domains. 



However, very recently, Benchimol et al. [4] suggested a binary heuristic, 
based on the MST relaxation, that we call RemoveMaxRC. It consists in remov- 
ing from Gp the tree arc of maximum replacement cost, i.e. the arc which would 
involve the highest cost augmentation if it was removed. By tree arc, we mean 
the fact that it appears in the MST of the last iteration of the Held and Karp 
procedure. Acutally, as shown in section [6j this branching leads to poor results 
and should not be used. 

Finally, Focacci et al. solve the TSPTW 15] by guiding the search with time 
windows, which means that the efficiency of CP for solving the ATSP should 
not rely entirely on its branching heuristic. 

4 Considering the reduced graph 

In this section, we consider a subproblem which is not a subset of constraints, 
as usual, but consists in the whole ATSP itself applied to a more restrictive 
scope: the reduced graph of Gp. The structure of the reduced graph has already 
been considered in a similar way for path partitioning problems |3l5j . In this sec- 
tion, we first study structural properties that arise from considering the reduced 
graph. Second, we show how to adapt such information to some state of the art 
implied models, including cost based reasonings. 

4.1 Structural properties 

We introduce a propagator, the Reduced Path propagator, which makes the 
reduced graph a (Hamiltonian) simple path and ensures by the way that each 
node is reachable from s and can reach e. It is a monotonic generalization of the 
algorithm depicted in [ST\ . Necessary conditions for this propagator have already 
been partially highlighted in [5], We first modify them in order to fit with the 
TSP and our model. Next, we provide a linear time incremental algorithm. 

Definition 1. Reduced path guarantees that any arc in Gp that connects two 
SCC, is part of a simple path which go through every SCC of Gp. 

Proposition 1. Given any directed graph G, its reduced graph Gp contains at 
most one Hamiltonian path. 

Proof. Let us consider a graph G such that Gp contains at least two Hamiltonian 
paths pi and pi, pi ^ p2. Since both pi and p2 are Hamiltonian then there exists 
at least two nodes {x, y} C Vp, x ^ y, such that x is visited before y in pi and 
after y in pi. Thus, the graph P = pl(Jp2 contains a path from x to y and 
from y to x. This is a circuit. As P C Gp, Gp also contains a circuit, which is 
a contradiction. □ 

We note Gp the reduced graph of Gp. We remark that, as s has no pre- 
decessor then its SCC is the node s itself. Also, as e has no successor then 
sccOf (e)= {e}. To distinguish nodes of V from nodes of the reduced graph, we 
note sccOf (s) = sp and sccOf (e)= ep. It follows that any Hamiltonian path in 
Gp will be from sp to ep. 



Proposition 2. If there exists a Hamiltonian path from s to e in Gp then there 
exists a Hamiltonian path in Gr. 



Proof. Lets suppose that Gr has no Hamiltonian path from sr to ep. Then for 
any path pp in Gr starting at sr and ending at Br, there exist at least one node 
x G Vr, which is not visited by pr. Thus, for any path pe in Gp starting at s 
and ending at e, there exist at least one SCC x G Vr which is not traversed by 
Pe, so Vm G nodesOf (x), then u ^ pe- Thus any path in Gp starting at s and 
ending at e is not Hamiltonian. □ 

It follows that any transitive arc of Gr must be pruned and that remaining 
arcs of Gr are mandatory (otherwise the graph becomes disconnected): any SCC, 
but Br, must have exactly one outgoing arc. An example is given in figure [T] 
the graph Gp contains four SCC. Its reduced graph, Gr, has a unique Hamilto- 
nian path P R = ({A}, {£?, C}, {E, D, F}, {G}). Arcs of G R \P R are infeasible so 
(A, E) and (C, G) must be pruned from Gp. 



Fig. 1. Reduced Path filtering, transitive arcs (dotted) are infeasible. 

We introduce a new data structure in Gr that we call outArcs : for each 
node x G Vr, outArcs(:r) is the list of arcs {{u,v) G A P \ sccOf(u)= x and 
sccOf (v)^ x}. We can now easily draw a complete filtering algorithm for the 
Reduced Path propagator which ensures the GAC over the property that Gr 
must be a path in 0(n + m) time: 

1. Data structures: Compute the SCC of Gp (with Tarjan's algorithm [32]) and 
build the reduced graph Gr = (Vr, A r ). 

2. Taking mandatory arcs into account: \/(u,v) G Am such that x =sccOf(«) 
and x ^sccOf(v), V(p, q) G outArcs(a;) \{(m,d)} remove arc (p,q). 




(a) G P 



(b) Gft before filtering 



(c) Gr after filtering 



3. Consistency checking: Make Gr a path if possible, fail otherwise. 

4. for each arc (u.v) G Ap such that x =sccOf (u) and y =sccOf (v), x =/= y, 

(a) Pruning: if (x,y) ^ Ar, remove arc (u, v). 

(b) Enforcing: if (x,y) € Ar and (u,v) is the only arc of Ap that links x 
and y, enforce arc (w,i>). 

A procedure performing step [3] starts on node sr, finds its right successor next 
(the one which has only one predecessor) and removes other outgoing arcs. Then, 
the same procedure is applied on next and so on, until ep is reached. Such an 
algorithm must be performed once during the initial propagation. Then, the 
propagator reacts to fine events. To have an incremental behavior, the propagator 
must maintain the set of SCC and the reduced graph. Haeupler et al. [H] worked 
on maintaining SCC and a topological ordering of nodes in the reduced graph, 
but under the addition of arcs. We deal with arc deletions. Moreover, we may 
have lots of arc deletions per propagation (at least for the first ones), thus we 
should not use a completely on-line algorithm. 

SCC maintenance: Let us consider an arc (u, v) £ Ap such that sccOf (u)= x 
and sccQf(?;)= y. If x ^ y and if (u, v) is removed from Ap then it must be 
removed from outArcs(:r) also. If x = y then the removal of (it, v) may split the 
SCC x, so computation is required. As many arcs can be removed from one prop- 
agation to another, we suggest a partially incremental algorithm which computes 
exclusively SCC that contains at least one removed arc, each one exactly once. 
We introduce a bit set to mark nodes of Gr. Initially, each node is unmarked, 
then when a removed arc, inside an unmarked SCC x, is considered, we apply 
Tarjan's algorithm on Gp H nodesOf(a;) and mark x. Tarjan's algorithm will 
return either x if x is still strongly connected, or the new set of SCC induced 
by all arc removals from x. In both cases, we can ignore other arcs that have 
been removed from nodesOf (x). Since the SCC of a graph are node-disjoint, the 
overall processing time of a propagation dealing with k arc deletions involving 
some SCC K C Vr is YlxeK 0( n x + m x) = 0(n + m). 

Gr maintenance and filtering: Algorithm [T] shows how to get an incremental 
propagator that reacts to SCC splits. When a SCC x is split into k new SCC, 
the reduced graph gets k — 1 new nodes and must be turned into a path while 
some filtering may be performed on GV. The good thing is that there is no need 
to consider the entire graph. We note X C Vr the set of nodes induced by the 
breaking of SCC x. Since Gr was a path at the end of the last propagation, 
we call p the predecessor of x in Gr and s its successor. Then, we only need to 
consider nodes of X \J{p 7 s} in Gr. To compute new arcs in Gr it is necessary to 
examine arcs of Gp, but only outArcs(p) and arcs that have the tail in a SCC 
of X need to be considered. Note that we filter during the maintenance process. 

Once we get those data structures, then it is worth exploiting them the most 
we can, to make such a computation profitable. Especially, a few trivial ad hoc 
rules come when considering SCC. We call an indoor a node with a predecessor 
outside its own SCC, an outdoor, a node with a successor outside its SCC and a 
door a node which is an indoor or/and an outdoor. 



Algorithm 1 Incremental Reduced Path Propagator 



Let x be the old SCC split into a set X of new SCC 

p <— GR.predecessors(x).f irst() {get the first (and unique) predecessor of x in Gr} 
s <— GR.successors(x).f irst() 
if (VISIT(p, s) / |X| + 2) then 

FAIL 
end if 



int VISIT(int current, int last) 
if (current — last) then 

return 1 
end if 
next 4 1 

for (node x £ current. successors) do 
if (|x.predecessors| — 1) then 
if (next ^ — l) then 

return {next and x arc incomparable which is a contradiction} 
end if 
next •<— x 
else 

GR.removeArc(current, x) 
end if 
end for 

for (arc (u. v) £ outArcs(current)) do 
if (sccOf (v) 7^ next) then 

Gp.removeArc(u, v) {Prune infeasible arcs} 
outArcs.remove(u, v) 
end if 
end for 

if (|outArcs(current)| — 1) then 

GM.addArc(outArcs(current).getFirst()) {Enforce mandatory arcs} 
end if 

return 1 + VISIT(next, last) 



Proposition 3. If a SCC X has only one indoor i G X, then any arc (j, i) G X 

is infeasible. 

Proof. First we remark that i cannot be s since s has no predecessors. Let us 
then suppose that such arc (j, i) G X is enforced. As the TSP requires nodes 
of to have exactly one predecessors, all other predecessors of i will be 

pruned. As i was the only indoor of X, then X is not reachable anymore from 
which is by Proposition [2] a contradiction. □ 

By symmetry, if a SCC X has only one outdoor i G X, then any arc (i, j) G X 
is infeasible. Moreover, if a SCC X of more than two nodes has only two doors 
i,j G X, then arcs and are infeasible. 

4.2 Strengthening other models 

In general, the reduced graph provides three kinds of information: Precedences 
between nodes of distinct SCC; Reachability between nodes of the graph; Car- 
dinality sets G VR\{eji}, |outArcs(:z;)| = 1. Such information can be directly 
used to generate lazy clauses [T3]. It can also improve the quality of the channel- 
ing between the graph variable and position variables by considering precedences: 
When adjusting bounds of position variables (or time windows), the BFS must 
be tuned accordingly, processing SCC one after the other. 



Some propagators such as DomReachability |26j . require the transitive clo- 
sure of the graph. Its computation requires 0{nm) worst case time in general, but 
since the reduced graph is now a path, we can sketch a trivial and optimal algo- 
rithm: For any node v £ V, we call S v C the set of nodes reachable from v 
in Gp and D v C Vr the set of nodes reachable from sccOf(v)<E Vr in Gr, includ- 
ing sccOf (v). More formally, S v = {u £ V\v — > u} and D v = xL){y £ Vr\x — > y}, 
where x =sccOf(u). Then, for any node v £ V, S v = {nodesOf (y) \y £ D v }\{v}. 
As Gr is a path, iterating on D v requires 0(|D. U |) operations. Also, since SCC 
are node-disjoints, computing S v takes 0(\D V \+J2 y( - D ^ |nodesOf (y)\) — 0(\S V \) 
because \D V \ < \S V \ + 1 and |{nodesOf (y) \y £ D v }\ = \S V \ + 1. As \S V \ < n, the 
computation of the transitive closure takes 0(^2 veV \S V \) which is bounded by 
0(n 2 ). It can be performed incrementally by considering SCC splits only. 

Finally we show how the MST relaxation of the TSP can be improved by 
considering the reduced graph. We call a Bounding Spanning Tree (BST) of Gp 
a spanning tree of Gp obtained by finding a minimum spanning tree in every 
SCC of Gp independently and then linking them together using the cheapest 
arcs: 

BST(Gp) = LU Gr MST(G P f| nodesOf (x)) 

Uaev R min f{( u > v ) \( u i v ) € outArcs(a)}. 

The resulting spanning tree provides a tighter bound than a MST. Indeed, 
since BST and MST both are spanning trees, f(BST(G P )) > f(MST(G P )), 
otherwise MST is not minimal. 

We will now see how to improve the Weighted Spanning Tree (WST) con- 
straint, leading to the Bounding Spanning Tree (BST) propagator. We assume 
that the reader is already familiar with this constraint, otherwise papers [29 30 
should be considered as references. The BST can replace the MST of the WST 
constraint: the pruning rules of WST constraint will provide more inference since 
the bound it tighter. Actually, we can do even better by slightly modifying the 
pruning rule of the WST constraint for arcs that are between two SCC: an 
arc linking two SCC can only replace (or be replaced by) another arc linking 
those two same SCC. Consider a BST of cost B, the upper bound of the ob- 
jective variable UB, an arc (x,y) £ Ar and a tree arc {u,v) £ outArcs(x), we 
can rephrase the pruning rule by: Any arc {u2,v 2 ) £ outArcs(a;) is infeasible if 
B — f(u, v) + f(u 2 , v 2 ) > UB. The reader should notice that no Lowest Common 
Ancestor (LCA) query is performed here. This do not only accelerate the algo- 
rithm, it also enables more pruning, because a LCA query could have returned 
an arc that does not link SCC x and y. Such an arc cannot replace (u, v) since 
exactly one arc of outArcs(a;) is mandatory. 

We now briefly describe a simple and efficient way to compute the BST. We 
assume that the Reduced Path propagator has been applied and that Gr is thus 
a path. Initially the BST is empty. First, we add to the BST mandatory arcs of 
Gp, then for each x £ Vr we add mirif((u,v) £ outArcs(x)). Finally, we run 
Kruskal's algorithm as described in |29|30j until the BST has n — 1 arcs. A faster 
way to compute a BST is to perform Prim's algorithm on successive SCC, but 
this method does not enable to use the efficient filtering algorithm of Regin [29 . 



Figure[2]illustrates this relaxation : the input directed graph, on figure 2(a) 



is 

composed of four SCC {A} 7 {B, C}, {E, D, F} and {G}. For simplicity purpose, 
costs are symmetric. Its minimum hamiltonian path, figure 2(b) costs 28 and 
we will suppose that such a value is the current upper bound of the objective 
variable. The MST of the graph, figure [2(c) only costs 19, which is unfortunately 
too low to filter any arc. Instead, the BST, figure 2(d) is much more accurate. 
It actually consists of the MST of each SCC, {0, {(BC)}, {(D, F), (E, F)}, 0} 
with respective costs {0, 10, 10, 0}, and the cheapest arcs that connect SCC each 
others: {(A, B), (C, D), (F,G)} with respective costs {2,3,2}. Thus, the entire 
BST costs 27. It is worth noticing that it enables to filter arcs (B, E) and (E, G). 
Indeed, (B, E) can only replace (C, D) in the relaxation, so its marginal cost is 
f(BST) + f(B, E) - f(C, D) = 27 + 5-3 = 29 which is strictly greater than 
the upper bound of the objective. The same reasoning enables to prune {E, G). 




5 The Held and Karp method 

The Lagrangian relaxation of Held and Karp has initially been defined for solving 
symmetric TSP. Instead of converting asymmetric instances into symmetric ones, 
through the transformation of Jonker and Volgenant [1] , we can directly adapt it 
to the asymmetric case: at each iteration fc, we define two penalties per node v € 
V, 7r-„ (v) and 7r^ t (^), respectively equal to (S'^ (v) - 1) * C (fe) and (<5+( fe ) (v) - 
1) * C {k \ We note the in-de gree of v in the MST of iteration k whereas 

5 + ( k \v) is its out-degree and is a constant whose calculation is discussed 
in [19|20| . As a path is expected, we post 7r|^(s) = 7Tott( e ) = 0- ^ rc costs are 
then changed according to: f^ k+1 \x,y) — f^ k \x,y) + ir^lx) + 7r-^(y). It has 
to be noticed that, since it relies on the computation of successive MST, such a 
model is equivalent to what would give the usual Held and Karp scheme used on 
a transformed instance. However, this framework is more general and can also 
handle the computation of Minimum Spanning Arborescence. 



Such a method should be implemented within a specific propagator to be 
easily plugged, or unplugged, into a constraint. We noticed that keeping track 
of Lagrangian multipliers from one propagation to another, even upon back- 
tracking, saves lots of computations and provides better results on average. Our 
approach is based on a few runs. A run is composed of K iterations in which we 
compute a MST according to Prim's algorithm and update and the cost 
matrix. Then, we run a Kruskal's based MST to apply the complete filtering 
of |4l29j . We first chose K = 0(n) but this led to disappointing results when 
scaling up to a hundred nodes. We thus decided to fix K to a constant. The value 
K = 30 appeared as a good compromise between speed and accuracy. Remark 
that, as we perform a fix point, the method may be called several times per 
search node, and since it is relatively slow, we always schedule this propagator 
at the end of the propagation queue. 

This procedure has the inconvenient of not being monotonic 1 (it is not even 
idempotent): filtering, related to other propagators, can slow down the conver- 
gence of the method. The intuition is that to go from a MST to an optimal tour, 
it may be easier to use some infeasible arcs during the convergence process. One 
can see an analogy with local search techniques that explore infeasible solu- 
tions in order to reach the best (feasible) ones more quickly [5]. This fact, which 
occured during some of our experiments involving static branching heuristics, 
breaks the usual saying the more filtering, the better. Moreover, it follows that 
we cannot measure precisely the improvement stemming from additional struc- 
tural filtering. We mention that the BST relaxation can be used within the Held 
and Karp scheme, however, this may also affect the convergence of the method 
and thus sometimes yield to poorer results. For that reason, we recommend to 
use a Lagrangian BST relaxation in addition to, rather than in replacement of, 
the usual Held and Karp procedure. 

6 Experimental study 

This section presents some experiments we made in order to measure the impact 
of the graph structure. We will show that branching according to graph structure 
only outperforms current state of the art results while using implied filtering 
based on graph structure avoids pathological behaviors on hard instances at a 
negligible time consumption. Our implementations have been done within the 
CHOCD solver which is an open source Java library. Tests have been performed on 
a Macbook pro under OS X 10.7.2 and with a 2.7 GHz Intel core i7 and 8Go of 
DDR3. We set a limit of 3 Go to be allocated to the JVM. We tested TSP and 
ATSP instance of the TSPLIB. For each one, we refer to the number of search 
nodes by |nodes| and report time measurements in seconds. As in [4], we study 
optimality proof and thus provide the optimal value as an upper bound of the 

1 A propagator P, involving a graph variable GV and a filtering function / : GV H > 
GV is said to be monotonic 31 iff for any GV' C GV,f(GV') C f(GV), where 
GV' C GV <s> G' P C Gp A G M C G' m . 



problem. We computed equivalent state of the art results (SOTA) (referred as 1- 
tree with filtering in [3]), to position our model in general. Their implementation 
is in C++ and has no memory restriction. 

Our implementation (referred as BASIC) involves one graph variable, one 
integer variable (the objective) and one single constraint that is composed of 
several propagators. Subtour elimination is performed by a special purpose in- 
cremental propagator, inspired from the NoCycle constraint [7]. The degree con- 
straint is ensured by special purpose incremental propagators described in section 
|3.1| The objective is adjusted by the natural relaxation and an implied prop- 
agator, based on the Held and Karp method. We mention that we solved rbg 
instances (that are highly asymmetric) by replacing the tree based relaxation 
by a Minimum Assignment Problem relaxation (also in SOTA). For that, we 
have implemented a simple Hungarian algorithm. Indeed, it always provided the 
optimal value as a lower bound at the root node. When a relaxation finds an 
optimal solution, this one can be directly enforced [3]. However, it could be in 
contradiction with side constraints. Thus, we unplugged such a greedy mode. 
The solver works under a trailing environment. 



6.1 Dedicated heuristics 



We experimentally compare the branching heuristics RemoveMaxRC and Sparse 
of section [3?3l We also introduce three variants of these methods: 

- EnforceMaxRC, consists in enforcing the tree arc of maximum replacement 
cost. It is the opposite of RemoveMaxRC. 

- RemoveMaxMC, consists in removing the non tree arc of maximum marginal 
cost, i.e. the arc which would involve the highest cost augmentation if it was 
enforced. This heuristic may require an important number of decisions to solve 
the problem. There are low probabilities to make wrong decisions, but if a mis- 
take has been performed early in the search tree, it might be disastrous for the 
resolution. 

- EnforceSparse, which first selects the set of nodes X with no successor in Gm 
and the smallest set of successors in Gp. Second, it finds the node x € X which 
maximize J2(x y)€A P I{( z i2/) ^ Ap\z e X}\. Then it fixes the successor of x by 
enforcing the arc (x, y) g Ap such that \{(z, y) £ Ap\z S X}\ is maximal. 

All branching heuristics are performed in a binary tree search. Remove- 
MaxRC, RemoveMaxMC and Sparse can be said to be reduction heuristics. They 
respectively involve a worst case depth for the search tree of 0(n 2 ), 0(n 2 ) and 
O(nlogn). In contrast, EnforceMaxRC and EnforceSparse perform assignments, 
leading to a 0(n) depth in the worst case. Assignment heuristics perform strong 
decisions that bring more structure in left branches of the search tree while it 
is the opposite for reduction branchings that restrict more right branches. An 
exception is Sparse which has a balanced impact on both branches. 
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Table 1. Search heuristics comparison on ATSP instances from the TSPLIB, with a 
time limit (TL) of 1,800 seconds and a memory limit (ML) of 3 Go. 
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Table 2. Implied filtering comparison with EnforceSparse heuristic. Hard instances 



ftv!70 and p43 are significantly improved by all filtering algorithms. 



6.2 Structural filtering 



We then study the benefit we could get from adding some implied filtering algo- 
rithms to the BASIC model with RemoveMaxMC and Sparse heuristics: 

- ARB: Arborescence and AntiArborescence propagators used together. 

- POS: The model based on nodes position, with an AllDif f erent constraint 
performing bound consistency. 

- AD: AllDif f erent propagator, adapted to graph variables, with GAC. 

- BST: Reduced Path propagator with a Lagrangian relaxation based on a BST, 
in addition to the usual Held-Karp scheme. 

- ALL: combine all above propagators. 



6.3 Results and analysis 

Table [T] provides our experimental results over the impact of branching strategies. 
RemoveMaxRC can be seen as our implementation of the SOTA model. The 
main differences between these two are the fact that we perform a fixpoint and 
implementation details of the Held and Karp scheme. As it can be seen, our 
Java implementation is faster and more robust (brl7, krol24p and rbg323). 
Results clearly show that the most recently used heuristic [3] is actually not very 
appropriate and that the more natural EnforceMaxRC is much more efficient. 
EnforceMaxRC is in general better than RemoveMaxRC, showing the limits 
of the first fail principle. The worst heuristic is clearly RemoveMaxMC while 
the best ones are EnforceMaxRC, Sparse and EnforceSparse. More precisely, 
graph based heuristics are the best choice for the f tv set of instances whereas 
EnforceMaxRC behaves better on rbg instances. This efficiency (not a single 
failure for instances rbg) is explained by the fact that EnforceMaxRC is driven 
by the MAP relaxation, which is extremely accurate here. In contrast, Sparse 
heuristic does too many decisions (because of the dichotomic branching), even 
if the number of wrong ones is negligible. Results show that assignment based 
branchings work better than reduction heuristics. The most robust branching 
strategy is EnforceSparse. Thus, graph structure can play a significant role during 
the search. However, side constraints of real world applications might require to 
use dedicated heuristics, so results should not entirely rely on the branching 
strategy. 

We now study the impact of implied filtering algorithms on robustness of 
the model. For that, we consider the best heuristic, EnforceSparse, and solve 
TSPLIB's instances within different model configurations. Results can be seen 
on Table [2j It can be seen that implied algorithms do not significantly increase 
performances on all instances, but it seriously improves the resolution of hard 
ones (ftvl70 and p43 are solved 5 times faster by ALL). Moreover, those extra 
algorithms are not significantly time expensive. The eventual loss is more due 
to model instability rather than filtering algorithms' computing time. Indeed, a 
stronger filtering sometimes yields to more failures because it affects both the 
branching heuristic and the Langrangian relaxation's convergence. In general, no 
implied propagator outperforms others, it depends on instances and branching- 
heuristics. The combination of them (ALL) is not always the best model but it 
provides robustness at a good trade off between filtering quality and resolution 
time. 

6.4 Consequences on symmetric instances 

In this section, we show the repercussion of our study on the (symmetric) Trav- 
eling Salesman Problem (TSP), which can be seen as the undirected variant of 
the ATSP. For that, we use an undirected model, as in 0]: Each node has now 
two neighbors. Previously mentioned implied structural filtering algorithms are 
defined for directed graphs, and thus cannot be used for solving the TSP. How- 
ever, our study about search heuristics can be extended to the symmetric case. 



It is worth noticing that the Sparse heuristic cannot be used here because the 
dichotomic exploration is not defined for set variables that must take two val- 
ues. We suggest to measure the impact of EnforceSparse heuristic that we have 
introduced and which appeared as the best choice for solving the ATSP. As can 
be seen in Table [3j the EnforceSparse heuristic can dramatically enhance perfor- 
mances on TSP instances. While the SOTA model fails to solve most instances 
within the time limit of 30 minutes, our approach solves all instances that have 
up to 150 nodes in less than one minute (kroB150 appart) and close half of the 
others that have up to 300 nodes. This seems to be the new limit of CP. Scaling 
further on simple instances would be easy: the number of iterations in the La- 
grangian relaxation should be decreased to get reasonable computing times, but 
solving bigger and still hard instances would require serious improvements. 
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Table 3. Impact of EnforceSparse heuristic on medium size TSP instances of the 
TSPLIB. Comparison with state of the art best CP results (column SOTA) [4] with a 
time limit (TL) of 1,800 seconds. 



7 Conclusion 

We have provided a short survey over solving the ATSP in CP and shown how 
general graph properties, standing from the consideration of the reduced graph, 
could improve existing models, such as the Minimum Spanning Tree relaxation. 
As future work, this could be extended to scheduling oriented TSP (TSPTW for 
instance) since Reduced Path finds some sets of precedences in linear time. 

We also provided some implementation guidelines to have efficient algo- 
rithms, including the Held and Karp procedure. We have shown that our model 
outperforms the current state of the art CP model for solving both TSP and 
ATSP, pushing further the limit of CP. Our experiments enable us to state that 
graph structure has a serious impact on resolution: not only cost matters. More 
precisely the EnforceSparse heuristic provides impressive results while implied 
structural filtering improves robustness for a negligible time consumption. 



We also pointed out the fact that non monotonicity of the Lagrangian relax- 
ation could make implied filtering decrease performances. An interesting future 
work would be to introduce some afterglow into the Held and Karp method: 
when the tree based relaxation is applied, it first performs a few iterations al- 
lowing, but penalizing, the use of arcs that have been removed since the last call 
of the constraint. This smoothing could make the convergence easier and thus 
lead to better results. 
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